HIGH ACCURACY METHOD FOR INTEGRAL EQUATIONS WITH 
DISCONTINUOUS KERNELS 
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Abstract. A new highly accurate numerical approximation scheme based on a Gauss type 
Clenshaw-Curtis Quadrature for Fredholm integral equations of the second kind 

r" 

x(t) + / k(t, s)x(s)ds = y(t) 



whose kernel k(t, s) is either discontinuous or not smooth along the main diagonal, is presented. This 
scheme is of spectral accuracy when k(t, s) is infinitely diffcrcntiable away from the diagonal t = s, 
and is also applicable when k(t, s) is singular along the boundary, and at isolated points on the 
main diagonal. The corresponding composite rule is described. Application to integro-differential 
Schroedinger equations with non-local potentials is given. 
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1. Introduction. Let the integral operator, 

[Kx){t)= ( k(t,s)x(s)ds, a<t<b, 

J a 

map Cy a b j, q > 1, into itself. In the present paper, we consider the numerical solution 
of the corresponding Fredholm integral equation of the second kind, 

(1.1) x(t) + [ k(t,s)x(s)ds = y(t), yeC q , a<t<b. 



When the kernel k(t, s) has a discontinuity either by itself or in its partial derivatives 
along the main diagonal t = s, one can not expect a high accuracy numerical ap- 
proximation based on Newton-Cotes or Gaussian Quadratures, ( see e.g. Figure 2 of 
Section 6 ), since, except for the Trapezium rule, the standard error bounds for these 
rules are not applicable. In the present paper we introduce a high accuracy discretiza- 
tion technique based on a Gauss type Clenshaw-Curtis quadrature for a certain class 
of such kernels which we call semismooth. 

Definition 1. A kernel k(t,s) is called (pi,p 2 ) — semismooth, if 



k(t,s) 



ki(t, s) if a < s < t 
k 2 (t,s) if t<s<b, 



where ki(t,s) 6 C^ b]x[a b] and k 2 {t, s) £ C^ b]x[a b] for some Pl ,p 2 > 1. 

Note that for our purpose each of the auxiliary kernels ki(t, s) and k 2 (t, s) must be 
defined in the whole square [a, b] x [a, b]. The convergence of our method is of 0(n 1 ~' r ), 
where r — mh\{pi,p 2 , q}. When r = oo, the convergence is superalgebraic, or spectral. 
Even for some singular kernels, when the obtained error estimates are not applicable, 
the method still shows good accuracy on numerical examples due to the clustering 
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of mesh points around singularities. The 2-step method of Deferred Approach to 
the Limit, based on the Trapezium rule, ( see e.g. Baker ||, p 363), works well for 
nonsingular kernels, but it is much more time consuming, for comparable accuracy, 
and is not applicable to kernels with singularities. 

In a way the present method is an extension of a previously examined situation, 
which occurs when a second order differential equation, such as the Schroedinger equa- 
tion with local potentials, is transformed into an integral equation, The kernel 
of this integral equation is generally semismooth, but with the additional structure 
of ki^(t,s) being of low rank. Such kernels we call semiseparable. In this case 
the algorithms developed in g, are adequate and give fast and accurate solution. 
However, for Schroedinger equations with non-local potentials, the corresponding ker- 
nels may be semismooth but not semiseparable, for which our present technique is 
perfectly well suited. This situation is examined in more detail in Section 7. 



In Section 2, we describe the discretization of equation (1.1) which is based on the 
Clenshaw-Curtis quadrature for a smooth kernel k(t,s). This discretization is differ- 
ent from the usual Gauss-Chebyshev quadrature, (e.g. Delves-Mohamed, Q), even for 
smooth kernels, and from that of Reichel jl5| , based on Chebyshev polynomial expan- 
sions. In Section 3 we consider semismooth kernels and show that the application of 
Clenshaw-Curtis quadrature results in a linear system of equations whose coefficient 
matrix is defined in terms of Schur, or componentwise, products of given matrices. 
The accuracy of approximation is determined by the smoothness of k\ and ki only, 
and is not affected by the discontinuity along the diagonal t = s. For smooth ker- 
nels this discretization is identical with the one described in Section 2. The proposed 
method works well also for the case when a semi-smooth kernel has singularities on the 
boundary of the square [a, b] x [a, 6], although the accuracy is not spectral anymore. 
The success of this method in this case is due to the clustering of Chebyshev support 
points near the boundaries, where the singularities occur. In Section 4 we describe 
the corresponding composite rule and in Section 5 we apply it to kernels which have 
finite number of singularities on the main diagonal t — s. In Section 6 we describe 
numerical experiments and comparisons with relevant existing methods for kernels 
with various discontinuities and singularities. In Section 7 we apply the developed 
technique to the solution of radial Schroedinger integro-differential equations with a 
non-local potential. 

2. Discretization of a Smooth Kernel. Let k(t, s) be differentiable in t and 
s. Assume that k(tk,s)x(s) as a function of s can be expanded in a finite set of 
polynomials , i.e., 

n 

(2.1) k(t k ,s)x{s) =J2 a kJ T j(s), -1<*<1, 

3=0 

where Tj(s) — cos{jarccos{s)) 1 j = 0, 1, ...,n, are the Chebyshev polynomials. 
Without any loss of generality we assume for now that a = — 1 and b = 1 in equation 



(|l_lj). Let 

/ r n+l 
k{t k ,s)x(a)da = J2b kj T j (t 

- 1 3=0 

Clenshaw and Curtis Q] showed that 

[frfcO, Hi, bkn+l] T = Si [dkO, a kl, — , Ctfcn 
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is the so called left spectral integration matrix. Here [v] denotes the transpose of 
the column vector v. Since Tj{\) = 1, for j = 0, 1, ...,n, it follows that 



/l n-t-l 
k(t k ,s)x(s)ds = J2 b *i 
- 1 j=0 



— [1; 1] [frfco, bki, ■■■,bkn+i] T — [1, 1]Sl [a k0 , a>ki> o, kn ] T • 
Let Tfc, k = 0, 1, n, denote the zeros of T n+ i, viz., 

(2fc + lW 

Tfc = COS 

so that 



2(n + l) 



(2fc + l)j7T 

rj-(rfc) = cos 2 - , k,j = 0,l,...,n. 



Substituting s = t/^/c = 0,1,..., n, into (2T), we obtain that 

k(t k ,T )x(T ) 

k(tk,Tl)x{T\) 



am 




Ctkl 








Oikn 





fc(tfe,T ra )x(r n ) 

where C _1 is an inverse of discrete cosine transformation matrix C whose elements 
are specified by 

Cjy =T j ( Tk ), k,j = 0,l,...,n. 

The matrix C has orthogonal columns, that is, C T C = diag(n, ^, j). Therefore, 
C _1 = diag(^, —, f^)C T . By choosing t k in (2.1) to be Chebyshev points and by 
substituting t = r k into (1.1), we get 



y(T k ) = x(r k ) + [1, l]SiC 1 diag(k(T k ,T ),k(T k ,Ti), fe(r fc ,r„)) 
Introducing [do, cr„] = [1, 1, lJS^C -1 we can write 

2/( T fc) = [c r o,cri,...,cr n ]dm3(fc(T fe ,T ),/c(r fc ,r 1 ),...,fc(T fc ,T n )) 



z(to) 
ac(n) 



z(r n ) 



z(to) 
x(ti) 



x(r n ) 
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or equivalently, 



y(r k ) = [fc(r fc , r ), k(r k , n), k(r k , T n )]diag(a , a x , <j n ) 



Therefore the discretization of the equation (1.1) for the case a 
as follows, 



a: (to) 
x(n) 

x{r n ) 
-1 and 6=1 is 



(2.2) 
where 



[I + KD ff ] x = y, 



K = (fc(T i ,r i ))^ =0! 

= diag(c7 , a\, a n ), 
x = [x(r ),x(ri), ...,a;(T n )] T , 

y = [y{ro),y{n), ...,y(r„)] T . 



The formulas (2.2) can be generalized for interval [a, b] other than [-1, 1] by the linear 
change of the variable h(r) = — a)r+ |(a+6). Thus if T)j=h(Tj), j = 0, 1, n, 
we have 



I+^KD. 



x = y, 



where 



K = (HVi,Vj))i,j=o^ D o- =diaflf(a-o,o-i,.„,cr„), 
x = [a:(??o),a;(7?i), x(?7„)] T , y = [y(r?o), y(Vi), ■■■,y(Vn)Y 



The accuracy of this discretization when k and x are not polynomials is discussed in 
a more general setting in the next section. 

3. Gauss Type Quadrature for a Semismooth Kernel. We consider now 
more general semismooth kernels, as in Definition 1, for which we write 

(3.1) x{t) + [ k 1 (t,s)x(s)ds+ [ k 2 (t,s)x(s)ds = y(t), a < t < b. 



In this section we describe the numerical technique for discretizing the equation (3.1). 
It is based on the Clenshaw-Curtis quadrature described in Section 2, which is well 
suited for computing antiderivatives. Let 

F(t)= J k 1 (t,s)x(s)ds, G(t,X)= J ki(t,s)x(s)ds, 

such that F(t) = G(t,t), and let 



m = f, 



k2(t, s)x(s)ds, J(t,\)= / k%(t, s)x(s)ds 



i 



A 
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Further, assume that fci(tjt, s)x(s) can be expanded in a finite set of polynomials, i.e. 
ki(tk, s)x(s) — Yl"— a ki Ti(s). As we have seen in Section 2, if 



(3.2) 

then 
(3.3) 



n+l 
3=0 



[PkOi fikl, Pkn+l] T — Sl [«feo, afel, afen] T • 



Similarly, assume that fe^fej s)x(s) = ^2™ =0 &kjTj(s). If 

~1 n+l 

J{t k ,X)= I k 2 (t k ,s)x{s)ds = y^^fcjTj(A), 

«* A J n 



then 



(3.4) 










• •,Pfcn+l 






IfifcO 


5*1, 
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is the right spectral integration matrix. Let r k ,k = 0, l,...,n, denote the zeros of 
T n+ \. Substituting A = Tfc, k — 0, 1, n, into (3.2), we obtain that 



G(tfc,7ft) 

G(t fc ,n) 

G(tfe,T„) 

and, similarly, 

J(tk,T ) 
J(tk,n) 

J(tk,T n ) 



CS L C 1 diag(k 1 (t k ,T ), kifa, r„)) 



CS_rC 1 diag(k 2 (t k ,T ), k 2 (t k ,T n )) 



. »(T n ) 

5C(T0) 
x(T n ) 



'n+l 



We remark that in writing the equality sign in (3.2) and (3.4), we assume that /?. 
is set to zero. This is an acceptable assumptio n b ecause G(t k ,X) is itself only ap- 
proximately represented by the polynomials in (|3.2| ) and the overall accuracy is not 
affected. Since F(r k ) = G(T k ,T k ) we get 



F{r k ) = [0, ...,0, 1,0, 0] CSLC-^diagihir^To), fci(r fe ,r„)) 



a (jo) 

z(T n ) 
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[wko, Wki, -, Wkn] diag(kt(T k ,T ), ki(r k , r„)) 



x(r n ) 



= [wko, Wki, —,Wkn\ diag(x(T )i ...,x{r n )) 



ki (T fc ,r ) 



ki (T fc ,r„) 



where [wko, Wkn] is the (fc + l)-st row of the matrix W ^ f CS^C -1 . We need now 
the following identity which can be verified by direct calculation. 

Lemma 3.1. Let A and B be n x n matrices and c — [ci, c n ] T . Then 
(A o B)c = diag(Adiag(ci, c„)B T ), where A o B denotes the Schur product of A 
and B, (A o B)^ = aijbij, i, j = 1, ...,n. 
Using this lemma we find that, 



(3.5) 



F(ro) 
Fin) 



diag(Wdiag(x(T ), ...,x(r n ))KT) = (WoKj) 

F(r n ) 

where Ki = (fci(ri, Tj))fj =0 . Similarly, 



x(t ) 
x(r n ) 



(3.6) 



H(r ) 
H(r n ) 



= (V o K 2 



x(r n ) 



where V = CS^C -1 . The formulas (3.5) and (3.6) can be generalized for an interval 
[a, b] other than [—1, 1] by the linear change of variables, h(r) — j(b — a)r + i(a + b). 
Thus if rjj — h{rj) 1 j = 0, 1, n, and with the notation 



F a (t) = / k 1 {t,s)x(s)ds, H b (t) 



we have 



(3.7) 



and, 



(3.8) 



F a (r)o) 

F a (m) 

F a {g n ) 
H b (r) ) 

H b (m) 



b — a 



-(WoKi) 



b — a 



(VoK 2 ) 



ki {t, s)x{s)ds, 



x(vo) 
x(vi) 

x(r]n) 

x(Vo) 
x(vi) 

X(f] n ) 



Using (3.7) and (3.8) we can now discretize the equation ( |3.lD as follows, 

b — a , 



(3.9) 



■(WoKi + VoK 2 ) 



x = y, 



7 



where x = [x(t]q), x(i] n )] T and y = [y(rjo), y(i] n )] T . Next we show th at i f 
the kernel function k(t, s) is smooth, such that ki = k 2 , then the discretization (3.9) 
reduces to (2.2). 

Proposition 1. Suppose that k(t,s) e C[ ab ] x [ ab ], an d that k\(t,s) = k 2 (t,s) = 
k(t,s). Then, 



I + ^£(WoKi + VoK 2 ; 



Proof. Without any loss of generality we assume that a = — 1 and 6 = 1. For 
t = t k , -l<t fc <l, 

(3.10) a;(f fc ) + G(t fc , A) + J(t k ,\) = y(t k ) 

holds for any A, — 1 < A < 1. Therefore if A = 1, then J(t k , 1) = and 

"fcO 

G{t k ,l)=J2Pkj =[1,...,1]S £ ' 



3=0 



OLkn 



[1, 1]S L C 1 diag(k(t k , r ), k(t k ,T n )) 



[o- ,ax, ...,a n ]diag(k(t k ,T ), ...,k(t k ,r n )) 



x(t ) 

x(r n ) 
s (t ) 

x(r n ) 



= [k(t k ,T ), k(t k , n), fc(ife, T„)]diag(cr , 01, cr n ) 

Substituting = r k for fc = 0, 1, ...,n, into ( 3.10| ), we obtain that 

[I + KDj] x = y. 

The assertion now follows. □ 



x(t ) 
x(n) 

x(r n ) 



We compared the numerical behavior of the discretization {2/2) and (3J3) for a 
number of smooth kernels and found that numerical answers differed in accuracy at 
the level of machine precision only. 

We now estimate the accur acy of approximation of the integral equation (3J) with 
the linear system of equations (3.E). The following property of Chebyshev expansions 
can be derived along the lines of an argument in Gottlieb and Orszag ( flQfl , p. 29). 

Proposition 2. Let f e C r [-1, 1], r > 1, and let 

oo 

/(t) = X>i2i(t), 1 < t < 1- 

3=0 



Then 



<*j\<- 

7T 



-i7rf{cos6) 
d9 r 
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and 



3=0 



< 



r — 1 n T 



It implies that if f(r) is analytic then the convergence of Chebyshev expansions 
is superalgebraic. Let now Fi(x) — f_-.f(t)dt and F r (x) — f f(t)dt. The 
following result can be found in Greengard and Rokhlin . 

PROPOSITION 3. Suppose that f £ CT_ X y, r > 1, and that f = (/(to), /(r„)) T , 
is i/ie vector of the function values at the roots of T n +i(x). Suppose further that Fi 
and F r are defined by 



Ft = (F l (T ),... 1 F l (T n )f 1 F r = (F r (r ),...,F r (r„)) 5 



The 



\Fi - CS L C- 



n r 1 



||F r -CS fl C~ 1 /|| 00 = 0(— ). 

n r 1 

Furthermore, all elements of the matrix CS^C -1 and CS/?C _1 are strictly positive. 



Let now rji = 



b-c 



where Tj is a zero of T n+ i(x), for i = 0, 1, 



,n, 

be the shifted Chebyshev points, and x = (x(r] ),x(rii), ...,x(r] n )) T be the vector 
of values of solution x(t) of equation (3.1) at rji. The following proposition follows 
immediately from standard properties of the Riemann integral, (see e.g. |T^] , p. 105). 

Proposition 4. Let k(t,s) be (pi,p 2 )-semismooth and let y(t) £ C? t ,, swc/i 
i/iai r = min{pi,p2, <?} > 1- Let the equation (3.1) define an invertible operator on 
C r [a b] . Then x £ C r [a b] . 

Let now 



and 



F, 



(Faivo), ■■■,F a (r] n )) 1 



H b = (H b (r) ),...,H b (r) n )) T . 
It follows from Proposition 3 that in conditions of Proposition 4, 

ll^-^(WoK 1 )x|| oc = 0(-L-), 



and 



U//,-^(VoK 2 )x|| oo = 0(^- T ). 



2 n- 
Combining the above results we obtain the following estimate for the residual. 



9 



Theorem 3.2. Let x be a solution vector of the equation i\3.{\ ), and x the vector 
of values of the solution x(t) att = r\j ll i — 0, 1, ...,n. Suppose that k{ t, s) is (pi,p 2 )- 
semismooth, and that y(t) £E C^ a b y Suppose further that the equation (3.1) defines an 
invertible operator on CT where r = min{pi , p2 , q} > 1. Then, 

||(I+^(WoK 1+ VoK 2 ))(x-x)|| 0O = 0(-l T ). 



It follows from the collectively compact operator theory, see Anselone [Q, that 
for sufficiently large n the matrices I + ^p(W o Ki + V o Ka), which depend on n, 
are invertible and their inverses are uniformly bounded. Therefore Theorem 1 implies 
that for increasing n, the convergence of a; to £ is of order 0(n 1 ~ r ). If p\ = p 2 = q = 
oo, then the convergence is superalgebraic. Numerical examples in Section 5 indeed 
demonstrate this type of convergence. 

4. The Composite Rule. In t his section, we describe the composite rule cor- 
responding to the quadrature of ( |3.9| ). Let 

a = bo < b\ < • • • < b m = b 

be a partition of the interval [a, b] , and let 



0,1, 



be the Chebyshev support points mapped into [bj—i,bj]. Define 
x(t)-- 



x x {t) 
x 2 (t) 



if b < t < b x 
if bi < t < b 2 



and 



y(t) 



(t) if b m -i <t< b m , 

Vi(t) if b <t<6i 
1ft (t) if h<t<b 2 

Vrn 

CO if &m-l <t<b m , 



and rewrite the equation (1.1) as a system of m equations, for j = 1, m, 



6o 



x j{t) + / fci(t, s)a!i(s)ds + • • • + / &i(t, s)xj(s)ds + / k 2 {t, s)Xj{s)ds 



bj-i 



(4.1) 



k 2 (t,s)x rn (s)ds = yj(t). 



Applying the quadrature of (3.S) to each of the integrals we obtain a system of linear 
equations as follows, for j = l,...,m, 



bi - b 



[(W + V)oKy]ii 



&i— i 



(WoK^VoKj 



&m-l 



[(W + V)oK mi ]4 
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where 







x{t[ 3) 




= (fci(^' ) , 




Kij 


= (fci(^'\ 






= (fc 2 (r«, 


'q 77 



))p,g=0. «/ 
))pfg=0> */ * >J. 



■,!/(r#)] T 



or in a block matrix form, 



(4.2) 



where 



A 2 i 



A12 
A 22 



^m2 



Ai 7n 

A2m 





Xi 








X2 















A^ = [I 



6,_i 



2 (WoK^lVoKj 



A - 6 J ~ 



[(W + V)oK Jt ], »/ i^j. 



In this paper we do not consider the issue of how to partition the interval [a, b] . An 
adaptive quadrature rule is possible here along the same lines as in Q, [ jl4| , namely, 
by using the size of last Chebyshev coefficients of k\ , k-2 and y in a given subinterval 
of partition to determine whether this subinterval should be further subdivided. This 
adaptive rule is a part of our research project in which we are going to compare the 
algorithm of the present paper with existing algorithms for Schroedinger equations 
with physically realistic nonlocal potentials. 

In general, the matrix ( |4.2| ) is not structured, and is being solved by the standard 
Gaussian elimination at the cost of 0(m 3 ) arithmetic operations, (we assume here that 
to is much larger than the nj's). If however, the semismooth kernel k(t, s) has some 
additional structure, then this structure is usually inherited by the matrix in (4.2). For 
example, if k\ and &2 are low rank kernels then the matrix A becomes semiseparable, 
and can be solved by existing linear complexity algorithms. We remark that in the 
case of the Schroedinger equation with non-local potentials discussed in Section 7 
below, the overall kernel is obtained as the composition of the semi-separable Green's 
function with the non-local potential. If the non-local potential is also semi-separable, 
which is the case when the non-locality arises from exchange terms, then the overall 
kernel is semi-separable as well, and the numerical techniques presented here, although 
still applicable, can be replaced by the IEM methods (|^], previously developed 
for local potentials. These methods give highly accurate linear complexity algorithms 
for the integral equation itself. 

If the kernel k(t, s) depends on the difference of the arguments, 



k(t,s) = k(\t - s\) = 



ki(t -s) if < s < t 



k 2 (t-s) if t < s < T, 
and if we use a uniform partition with the same number of points per partition, then, 
fcr(r»,rf )=k r (r p -r q ), r = l,2, 



11 



A: 


A 2 


A 3 • 


A 


A 2 


Ai 


A 2 • 


A m _i 


A 3 


A 2 


Ai • 


• A m _ 2 



Xi 




yi 


X 2 




m 









and we obtain a block Toeplitz matrix 



(4.3) 



A m A m -1 • ' ' A 2 Ai 

which can be solved in 0(m log(m)) or 0(m . ) arithmetic operations by algorithms 
available in the literature. 

Finally we would like to point out that the above composite rule can be used to 
handle kernels which have a hnite number of singularities on the main diagonal, t = s. 
In this case one has to include all the singular points as a subset of all partition points. 
This is illustrated in the next section on a simplest case of one such singularity, but 
in full detail. 

5. Kernels with Singularities on the Main Diagonal. Suppose that the 
kernel kit, s) has a singularity at (c, c), a < c < b. inside the square [a, b] x [a, 6], Since 
the Chebyshev points % = + s ^-, i = 0, 1,2, ...,n, are clustered towards the 

boundaries of the square [a, b] x [a, b] , we do not have sufficient values of k(t, s) around 
singular point (c, c) for a good approximation. Therefore we apply the composite rule 
of Section 4 with c being a partition point. For the sake of simplicity we assume that 
c is the only partition point, thus [a, b) is partitioned into two subintervals [a, c] and 
[c, b) . Without loss of generality, we consider the solution of 



(5.1) 
where 



x(t) 



k(t,s) 



k(t, s)x(s)ds = y(t), 



h(t,s) if -l<s<f 
k 2 (t, s) if t < s < 1, 



and assume that the kernel k(t, s) has a singular point at (0,0). Define 

T(ri f xi(t) if -l<t<0 
[ ' ~ \ x 2 (t) if < t < 1 



and 



y(t) 



2/i (t) if -1<*<0 
y 2 (t) if < t < 1. 



We can rewrite the equation (5.1) as a system of two equations. For — 1 < t < 0, 
(5.2) xi(t) + / k(t, s)xi(s)ds + / k(t, s)xi(s)ds + / k(t, s)x2(s)ds = yi(t), 



and for < t < 1, 



(5.3) x 2 (t) + / kit, s)xx{s)ds + / k(t, s)x 2 (s)ds + / k{t,s)x 2 (s)ds = y 2 {t). 
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Discretizations of the equations (5.2) and ( p.3[ ), respectively, are as follows. 

(5.4) [I + ~(W o K u + V o Ku)]x! + ~[(W + V) o K 2 i]5; 2 = y x 
and 

(5.5) ~[(W + V) o Kia]ii + [I + i(W o K 22 + V o K 22 )]x 2 = y 2 , 
where 

Ku = Mf'.f');,^. K n = fc 2 (r«,r«)^ =0 , K 12 = h(r^ , r^)^ 
K 22 = fcr(rf ,rf )^ =0 , K 22 = fc(r«,7f )£=<>, K 21 = fe^.f)^, 
with 



.(i) _ 



1 



.(2) _ 



i, i = 0,1, ...,n x , j = 0,1, ...,n 2 . 



-Ti - - and ' = -Tj 



Herexi = [^(rf ), ^(rV )] and x 2 = [a; 2 (rj; ),..., x 2 {t, 12 ')] 1 . In the matrix 
form we write, 



(5.6) 
where 





A12 " 










_ A 2i 


A 22 _ 




. * 2 







A n = 


[I + i(Wo 


K n + 


V 


oK n )] 


A i2 = 


i[(w + v) 


oK 21 






A 2i = 


i[(W + V) 


oK 12 






A 22 = 


[I + i(Wo 


K 22 + 


V 


oK 22 )] 



The size o f the matrix in the equation (5.6) is (ni+n 2 + 2) x (ni+n 2 + 2). The formulas 
(5.2) and (5^3) can be generalized for interval [a, c] and [c, b] other than [—1,1] by the 
linear change of the variable p(t) = ^(c — a)t+^(c + a) and q(t) = \{b — c)t + ^(b + c) 
if (c, c) is a singular point of k(t, s). Thus if 



q(n), i = 0,1,2, ...,m 



then. 



where 



" An 


A 12 




Xl 




yi 


A 2i 


A 22 




x 2 




y~2 



c — a 



b-c 
2 

c — a 



k-21 



(WoKn + VoKi 
z 

[(W + V)oK 21 ], 
[(W + V)oK 12 ], 



A 22 = [I + o K 22 + V o K 22 )], 

with Kij defined as above. A numerical example illustrating this technique is given 
in the next section. 
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6. Numerical Examples. In this section we compare our methods with some 
exiting algorithms for the following types of kernels, 

• Type 1 : Discontinuity along the diagonal t = s. 

Discontinuity in the first order partial derivatives along the diagonal 



Type 2 : 
t = s. 
Type 3 
Type 4 



Singularity on the boundary of the square and Type 2. 
Singularity on the main diagonal. 
These are the methods which have been implemented for comparison purposes, 
G-Leg : Nystrom discretization based on Gauss-Legendre quadrature. 

Two step Deferred Approach to the limit. Approximate solutions Xi, X2, 
and x% for subintervals of partition h, ~, and j, respectively, are computed. 
Then the numerical solution x(s) is obtained by ( see e.g. Baker pi) 



T-Def 



x{s) 



64x 3 (s) +xi(s) -20a; 2 (s) 
45 



Atk-T : Atkinson's iteration with the composite Trapezium rule. Applied to kernels 
k(t, s) which have discontinuities in the first order partial derivatives along 
the diagonal t = s. 



Alg-1 : Algorithm of Section 2, (2.2). 
Schur : Algorithm of Section 3, (3.E). 
Sch-C : Algorithm of Section 4, (O) 



The number of points used in discretizations is denoted by n. Error denotes 
| \x— a; T ||/||x||, where x and x T are the analytic and the numerical soluions, respectively. 
In each plot, log(Error) is the common logarithm of the Error. All computations were 
done on a DELL Workstation with operating system RedHat Linux 5.2 in double 
precision. All examples are set-up by choosing a simple analytic solution and then 
computing the corresponding right hand side. We remark that the values of x(t) are 
found inside the interval (or each of the subintervals of partition ) at Chebyshev points 
to,ti, ...,t„. The value of x{i) for t ^ T& can be found as follows. Applying C _1 we 
can find "Chebyshev- Fourier" coeffcients of x(t), 



a 




x(r ) 






x(n) 




= C 1 






_ x(r n ) 



Thus, 

n 

x(t) ^^ajTjihit)), a<t<b. 

3=0 

The value of Tj (t) for t ^ Tk is found now using the recursion satisfied by Chebyshev 
polynomials, T J+1 (t) = 2tT 3 {t) - T^i(t). 
Example 1. 

x(t) + \[ k(t,s)x(s)ds = y(t), -1 < t < 1, 



where y(t) = A(e + e" 1 ) + (1 - 2A)e _t , and 



k(t,s) = 



if 
if 



- 1 < s <t 
t < s < 1. 
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The analytical solution is x(t) = e _t . Since this kernel is discontinuous along the 
diagonal t = s, Gauss-Legendre quadrature gives low accuracy. The accuracy in the 
Atkinson's iteration improves very slowly. The algorithm of Section 3 gives accuracy 
of order 10 -15 with only 16 support points, whereas the 2-step method of Deferred 
Approach to the Limit method requires n = 256 points to achieve comparable accu- 
racy. Moreover, it requires computation of X2{t) and Xz{t) at the cost of 0((2n) 3 ) 
and 0((4n) 3 ), respectively. 

Fig 1 



-2 

-4 

-O -6 
o 

!±L -8 
o> 
o 

-10 
-12 
-14 



100 200 300 400 500 600 
Number of points 

Fig. 6.1. Comparison of numerical solutions of Example 1, A = 0.1. 



Example 2. 

x(t) + \[ k(t,s)x(s)ds = y(t), 0<t<T 
Jo 

where y{t) = (1 — Asm 2 ^ + \)sint + — t — sm ^ 2T ) ) \ eos ^ ailc j 




, . . ,. h / sm(t - s) if < s < t 
k(t, s) — sm( \t — s = < . ; J. .„ , . _ 
v ' ; Vl \ sm(s-t) if t < s < ~ . 

The analytical solution is x(t) = sin(t). This kernel has discontinuities in the first 
order partial derivatives along the diagonal, t — s. Again standard Nystrdm-type 
discretization methods fail to give high accuracy in this case. In the first experiment 
we take T = ^ and A = — — . Our method shows the order of 10~ 14 accuracy with 
only 16 points in [0, without any partitioning. The 2-step method of Deferred 
Approach to the Limit gives the accuracy of O(10~ 14 ) with n = 256, but at much 
higher cost than our method, see Fig 2. The second part of Example 2 is to compare 



the composite rule described in Section 5 with the basic quadrature (3.9) of Section 4 
when the length of the interval of integration, [a, 6], becomes increasingly large. Here 
M denotes the number of subintervals in [0, T] and Mn stands for the total number 
of support points in the given interval [0, T]. 
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Fig 2 






100 200 300 400 500 600 
Number of points 



Fig. 6.2. Comparison of numerical solutions of Example 2. 



Table 1 ( T = 200tt, A = -f ) 



M 


1 


1 


1 


2 


4 


8 


Mn 


256 


512 


1024 


256 


512 


1024 


Error 


2.4e + 01 


3.0e- 02 


CPUtimc limit exceed 


1.2e + 01 


7.8e - 02 


2.2e- 11 



Without partitioning, i.e. with M = 1, we increase the number of support points 
from n = 128 to n = 1024. For n = 512 the accuracy is of order O(10~ 2 ), but for 
n = 1024 the CPU time limit is exceeded. When the interval is partitioned into 8 
subintervals and n = 128 i.e., the total number of points is 1024, the accuracy now is 
of order <3(10- n ). 
Example 3. 

x(t) + J k(t,s)x(s)ds = y(t), -1 < t < 1 
where y(t) = 1 — t 2 + j^(a,rcta,n(t) — arctan(— 1)) — (1+t J 1+t ^^ , and 



[ (l_ t i)(l_ s 2) it t < s s i. 

The analytical solution is x(t) = 1 — t 2 . Since this kernel has singularities along the 
boundaries of the square [— 1, 1] x [— 1, 1] methods based on the Trapezium rule are not 
applicable. Therefore we compare our algorithms of Section 1 and Section 3 with the 
A^ysirom-Gauss-Legendre discretization only. The algorithm of Section 1 shows the 
same accuracy of numerical solution as the Gauss-Legendre quadrature. The method 
of Section 3 gives O(10~ 13 ) accuracy with n = 32 points, whereas Nystrom-GsMSS- 
Legendre quadrature gives O(10~ 3 ) with n = 256 points. 
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Fig 3 










-2 




_4 


o 


-6 


a_ 




L±J_ 




O) 
O 


-8 




-10 




-12 




-14 




50 



100 150 200 
Number of points 



250 



300 



Fig. 6.3. Comparison of numerical solutions of Example 3. 



Example 4. 



x{t)+ k(t,s)x(s)ds =y(t), -1 < t < 1 



where y(t) = 2(1 - t 2 + 2t 3 ) + (1 + 2t 4 )ln(t 2 + t 4 ) - ln(l + t 2 ) - 2tHn{\ + t 4 ), and 



k(t, s) 



l/(t 2 
l/(s 2 + t i ) 



s 4 ) 



if -l<s<t 
if t < s < 1. 



The analytical solution is x(t) — 4i 3 . The kernel k(t,s) has a singularity at (0,0). 



'(2i-l) 



Also is singular at t = 0. Since Chebyshev points cos( 2Af '' ttJ, « = 1, 2, n, are 
clustered towards the end points of interval [—1,1], discretization formula ( |3.9[ ) does 
not contain sufficient values of kernel near (t, s) = (0, 0). Therefore we partition [—1,1] 
into [—1,0] and [0, 1]. The choice of n — 256 Chebyshev points in each subinterval with 
the total of n — 512 points gives O(10 -11 ) accuracy. For comparison the best accuracy 
of the Gauss-Legendre quadrature without partitions and with n = 512 support points 
is O(10 -4 ) while the best accuracy of the algorithm of Section 2 without partitions 
is of O(10~ 6 ), see Fig 4. 



7. Application to Non-Local Schroedinger Equations. In this section we 
demonstrate that the developed numerical technique is also applicable to problems 
other than integral equations, for example, to integro-differential equations. We chose 
here the radial Schroedinger equation which models the quantum mechanical interac- 
tion between particles represented by spherically symmetric potentials. These poten- 
tials are usually local, i.e., they depend only on the distance between the two particles, 
in which case the equation is a differential equation which is routinely solved in com- 
putational physics. However, if there are more than two particles present, then the 
potentials can become non-local and the differential Schroedinger equation becomes 
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Fig 4 




100 200 300 400 500 600 
Number of points 



Fig. 6.4. Comparison of numerical solutions of Example 4- 



an integro-differential equation for the wave function ip, 

(7.1) + n 2 i,{r) = ( T v(r, r')i){r')dr' , 

dr 1 J 

which is defined for < r < oc, satisfies the condition ip(0) = 0, and is bounded 
at infinity. It is assumed that v(r,r') is negligible for r > T or r' > T, see e.g. 
[fl3[ . Because it is numerically more difficult to solve the Schroedinger equation in the 
presence of a nonlocal potential, the latter is customarily replaced by an approximate 
local equivalent potential. There is, however, a renewed interest in the nonlocal 
equations, and a significant number of papers on this subject appeared in the past 
few years, (our database search returned over 50 related publications). 



Using the technique of |8| it easy to show that ( 7_A ) is equivalent to the following 
integral equation, 

ip(r) H / sin(nr') / v(r' ,p)tlj(p)dpdr' 



u 



sin(Kr) 



T pT 

cos(nr') / v(r' ,p)ip(p)dpdr' = sin(Kr). 



(7.2) ^(r) + C ° s(Kr) [ T h(r,r')t/j(r')dr' + SlU ^ C k 2 (r,r')ip(r')dr' 



where 



sin(Kr), 



fci(r,r')= I sin(Kp)v(p 7 r')dp, k%{r,r') = f cos(np)v(p,r')dp. 
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We consider now the case when v(p, r') is (pi,p 2 )-semismooth, such that 

lP ' > ~ \ v 2 {p,r') if r' <p<T. 

In order to use the method which we developed in previous sections, we rewrite 
equation (^) as follows, 

i>{r) + — r ki(r,r')i>(r')dr' + ^ f k^r.r'^r^dr 1 

K JO K Jr 

(7.3) + T k 2 {r, r')i>(r')dr' + ^ f fc 2 (r, r')^{r')dr' = s(r), 



where for notational convenience we abbreviate, c(r) = cos(nr), s(r) = sin(nr). We 
have 

, J fen (r,r') if 0<r' <r 

Kl ^ r >- \ fc 12 ( r ,r') if 0<r<r', 



and 



fe 2 i(r,r') if r' <r <T 
k 22 (r,r') if r < r' < T, 



fe 2 (r,r') = 

where 

ku(r,r')= I s(p)vi(p,r')dp + J s(p)v 2 (p,r')dp, 

s(p)v 1 (p,r')dp+ / s(p)v 2 (p,r')dp - I s(p)v 2 (p,r')dp, 



k\ 2 {r,r')= j s(p)vi(p,r')dp 
Jo 

k 2 i(r,r')= c(p)v 2 (p,r')dp 

J r 

k 22 (r,r')= I c(p)v 1 (p,r')dp + J c(p)v 2 (p,r')dp, 

c(p)vi(p,r')dp- / c(p)v 1 (p,r)dp+ I c(p)v 2 (p,r')dp. 



Thus, 



^( r ) + £W f k n {r,r')ip(r')dr' + ^ / fei 2 (r, r')^(r')dr' 



(7.4) +^ I" k 21 (r,r'W(r')dr' + ^ I k 22 (r,r')^{r')dr' = s(r), 



Applying our quadrature to this equation we get, 

(7.5) [ I + ^D C (W o Ku + V o K 12 ) + ^B S (W o K 21 + V o K 22 ) $ = s, 
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where in more detail, 



1> 
D c 
D s 
D ff 
s 
W 
Kn 

K 12 

K21 



■ diag(cos(Kto), cos(nti), cos(Kt n )), 
- diag(sin(Kto), sm(reii), sin(Kt n )), 
= dia 5 ([l s ...,l]SiC- 1 ) ) 
= [sin(Kto), sin(Kti), sin(Kt n )] T , 

-cs L c-\ v = cs i? cr 1 , 

: (M*<.*j))«=o 

: |[diag(WD s (V! - V 2 )) + (WD S V 2 )] 

: (fci2(tj,tj))"j =0 

: |(WD a Vl), 



(fc2l(*i>*i))"j=0) 

= |(VD C V 2 ) 

K 22 = (k22(U,tj))™ j=0 

= |[(VD c Vi) + rf W5 (VD c (V 2 - Vi))]. 



We illustrate now the obtained discretization with two examples. In the first example 
we use a prototype of the Yukawa potential, (e.g. Jl^l , 23. c), which is simplified to a 
degree such that an analytic solution can be found. In our terminology this potential 
is semiseparable. We note once more that the case of this semi-separable potential 
could be treated more easily by the techniques already presented in || , and we use it 
here only because the comparison with the analytic solution is possible. 
Example 1. Let 



, ,, _ f \e p ~ r if < p < r' 
v[p,r)- I Xpr ,_ p . f r , < p < T 

It is easy to see that if i/>(r) = e~ r , then the right-hand side has the form, 



3Xk 

V{r) = (1 - — )e 



3Ak Xk 
— cos{r) - —re 



By comparing the analytical solution given above with the numerical solution of (7.5) 
at the discretization points, we get the following relative errors in the case of A = 
0.1, k = 1 and T = 20. 



n 


16 


32 


64 


128 


256 


Error 


1.2e+01 


3.4e - 07 


8.1e-09 


3.4e-09 


6.0e - 09 



In the second example we consider a more interesting case for which the techniques 
of H are not applicable. This time the non-locality is a prototype of the optical model 
Perey-Buck potential, (e.g. ||, Ch. V.2). In our terminology this potential is semi- 
smooth, but not semiseparable. 
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Example 2. Let 



v(p, r') = 




if < p < r' 
if r' <p<T. 



Solving ( [7.5| ) at n shifted Chebyshev support points t| , i = l,...,n, and 2n points 
sj 2 "\i = 1, ...,2n, we obtain numerical solutions ip^(r) and ?/>( 2n )(r), respectively. 

To get the values of ?/>( 2n )(r) at t\ n \ we follow the procedure described in the 
beginning of Section 6. The error e n is obtained by comparison of the solutions ij)^ n > 
and ■)/)( 2 ™) as follows, 

en = ||V (2n) (ti n) ) -^(t^JIU/H^^JHoo. 
Here we take A = 0.1, k = 1, A = 100, and T = 20. 



n 


8 


16 


32 


64 


128 


256 




l.Qe-0 


1.2e-03 


1.6e - 09 


7.7e - 15 


1.6e - 14 


4.8e- 14 



We see that for this choice of A the matrix is well conditioned and the double precision 
accuracy is obtained with 64 points. 

8. Summary and Conclusions. In this paper, which is one of a sequence treat- 
ing integral equations, we describe a new method for solving integral equations for 
the case when kernel can be discontinuous along the main diagonal. It has the fol- 
lowing advantages for a large class of such kernels: (i) for semismooth kernels it gives 
a much higher accuracy than it was ever possible with standard Gauss type quadra- 
ture rules; (ii) it is of comparable accuracy with Gauss type quadratures for smooth 
kernels; (iii) it exploits additional structure of the kernel such as a low semi-rank, or 
a displacement structure, for example, to allow for reduced complexity algorithms for 
the discretized equations, (iv) the numerical examples provided in the present study 
illustrate increased accuracy of our method compared to other more conventional 
methods. 

Our method is also applicable to other problems, such as the computation of 
eigenvalues and eigenfunctions of integral and differential operators and solution of 
integro-differential equations. 

Our method may find applications in quantum mechanical atomic and nuclear 
physics problems, where the requirement of indistinguishability of the electrons leads 
to non localities in the potential contained in the Schroedinger equation due to the 
presence of exchange terms. These, in turn, lead to integro-differential equations 
which are usually solved by iterative finite difference methods, or by orthogonal func- 
tion expansion methods. We plan to compare our new method with some of the 
existing methods in future investigations on more realistic examples. 
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